Configuraciones generales
# Paquetes
library(tidyverse)
library(lubridate)
library(readxl)
library(janitor)
library(mgcv)
library(gratia)
library(ggrepel)
set.seed(1)
dir.create("figs", showWarnings = FALSE)
# =================== Estética global ===================
theme_Ara <- function(base_size = 12, base_family = "") {
theme_minimal(base_size = base_size, base_family = base_family) %+replace%
theme(
plot.title = element_text(face = "bold", size = base_size + 2, hjust = 0, margin = margin(b = 6)),
plot.subtitle = element_text(color = "grey25", margin = margin(b = 8)),
plot.caption = element_text(color = "grey40", size = base_size - 2, hjust = 0),
axis.title = element_text(face = "bold"),
axis.text = element_text(color = "grey20"),
panel.grid.minor= element_blank(),
panel.grid.major= element_line(color = "grey88", linewidth = 0.4),
strip.text = element_text(face = "bold"),
legend.position = "top",
legend.title = element_text(face = "bold")
)
}
pal_okabe <- c("#0033A0", "#017574", "#54a2a4", "#D55E00", "#CC79A7", "#f6b617", "#00a2af", "#000000","#d92447")
ggsave_Ara <- function(filename, plot, width = 9, height = 5, dpi = 320) {
ggsave(filename, plot = plot, width = width, height = height, dpi = dpi, bg = "white")
}
# Predicciones con IC95% en escala de respuesta
pred_ci <- function(mod, newdata) {
p <- predict(mod, newdata = newdata, type = "link", se.fit = TRUE)
inv <- mod$family$linkinv
mu <- inv(p$fit)
lwr <- inv(p$fit - 1.96*p$se.fit)
upr <- inv(p$fit + 1.96*p$se.fit)
bind_cols(newdata, tibble(mu = mu, lwr = lwr, upr = upr))
}
# Curva de un smooth (usando gratia::smooth_estimates) + IC puntuales
smooth_curve <- function(mod, smooth,
x_candidates = c("x","t_meses","mes_num",".x"),
out_x = "x",
level = 0.95, linkinv = NULL) {
s <- gratia::smooth_estimates(mod, smooth = smooth)
est_col <- intersect(c("estimate",".estimate"), names(s))[1]
se_col <- intersect(c("se",".se","std.error",".std.error"), names(s))[1]
x_col <- intersect(x_candidates, names(s))[1]
stopifnot(!is.na(est_col), !is.na(se_col), !is.na(x_col))
if (is.null(linkinv)) linkinv <- mod$family$linkinv
z <- qnorm(0.5 + level/2)
s %>%
transmute(
!!out_x := .data[[x_col]],
est = .data[[est_col]],
se = .data[[se_col]],
lower = est - z*se,
upper = est + z*se,
mu = linkinv(est),
lwr = linkinv(lower),
upr = linkinv(upper)
)
}
# Segmentador de tramos contiguos a partir de bandera booleana
segmentar_tramos <- function(df, flag_col, gap_max = 1.5) {
df %>%
filter(.data[[flag_col]]) %>%
arrange(t_meses) %>%
mutate(gap = t_meses - lag(t_meses),
new_block = is.na(gap) | gap > gap_max,
tramo = cumsum(replace_na(new_block, TRUE))) %>%
group_by(tramo) %>%
summarise(
inicio = first(mes),
fin = last(mes),
n_meses = n(),
d_med = mean(derivative),
.groups = "drop"
)
}
theme_set(theme_Ara())
# ----- Switch para evaluar lluvia con AR(1) -----
usar_lluvia_ar <- TRUE # <— pon FALSE si no quieres correr la sección con lluvia
# Parámetros de entrada
ruta_xlsx <- "~/Library/Mobile Documents/com~apple~CloudDocs/CCGS/Manuscritos/Guacamayas_Lacandona/Ara_Lacandona_git/Bases_datos/BASE COMPLETA JUL2023_EMM.xlsx"
# Usa "BASE FINAL" si existe; si no, la primera hoja
sheets <- readxl::excel_sheets(ruta_xlsx)
sheet_use <- if ("BASE FINAL" %in% sheets) "BASE FINAL" else sheets[1]
raw <- read_excel(ruta_xlsx, sheet = sheet_use) |> clean_names()
# Columnas esperadas tras clean_names():
# fecha, transecto, no_de_transecto, tamano_de_grupo, no_adultos, no_juveniles, no_indeterminados, lluvia
dat <- raw |>
mutate(
fecha = suppressWarnings(as_date(fecha)),
mes = floor_date(fecha, "month"),
transecto = as.character(transecto) |> stringr::str_squish(),
no_transecto = suppressWarnings(as.integer(no_transecto)),
tam_grupo = suppressWarnings(as.numeric(tamano_de_grupo)),
n_adultos = suppressWarnings(as.numeric(no_adultos)),
n_juveniles = suppressWarnings(as.numeric(no_juveniles)),
n_indet = suppressWarnings(as.numeric(no_indeterminados)),
lluvia = suppressWarnings(as.numeric(lluvia))
) |>
select(mes, fecha, transecto, no_transecto,
lluvia, tam_grupo, n_adultos, n_juveniles, n_indet) |>
filter(!is.na(mes), !is.na(transecto))
# Agregación mensual por transecto
mensual_tr <- dat |>
group_by(mes, transecto, no_transecto) |>
summarise(
conteo = sum(tam_grupo, na.rm = TRUE),
adultos = sum(n_adultos, na.rm = TRUE),
juveniles = sum(n_juveniles, na.rm = TRUE),
indet = sum(n_indet, na.rm = TRUE),
lluvia = mean(lluvia, na.rm = TRUE),
.groups = "drop"
) |>
mutate(
mes_num = month(mes),
t_meses = as.numeric(difftime(mes, min(mes, na.rm = TRUE), units = "days"))/30.4375
)
# Total mensual
mensual_total <- mensual_tr |>
group_by(mes) |>
summarise(
conteo_total = sum(conteo, na.rm = TRUE),
adultos_total = sum(adultos, na.rm = TRUE),
juveniles_total = sum(juveniles, na.rm = TRUE),
lluvia = mean(lluvia, na.rm = TRUE),
.groups = "drop"
) |>
mutate(
mes_num = month(mes),
t_meses = as.numeric(difftime(mes, min(mes, na.rm = TRUE), units = "days"))/30.4375
)
# Vista rápida
dplyr::glimpse(mensual_total)
## Rows: 91
## Columns: 7
## $ mes <date> 2013-08-01, 2013-09-01, 2013-10-01, 2013-11-01, 2013-…
## $ conteo_total <dbl> 88, 128, 79, 89, 82, 114, 93, 63, 69, 130, 135, 116, 1…
## $ adultos_total <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ juveniles_total <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ lluvia <dbl> 535.70, 606.50, 647.60, 203.60, 311.00, 45.20, 56.00, …
## $ mes_num <dbl> 8, 9, 10, 11, 12, 1, 2, 3, 7, 8, 9, 10, 12, 1, 3, 6, 8…
## $ t_meses <dbl> 0.000000, 1.018480, 2.004107, 3.022587, 4.008214, 5.02…
# 2.1 GAM sin AR(1) para estimar rho inicial
gam_total0 <- gam(
conteo_total ~ s(t_meses, k = 20) + s(mes_num, bs = "cc", k = 12),
data = mensual_total,
family = nb(), method = "REML",
knots = list(mes_num = c(0.5, 12.5)),
na.action = na.exclude
)
# 2.2 Estimar rho de ACF(1) a partir de residuales (aprox.)
acf1 <- try(acf(residuals(gam_total0, type = "pearson"),
na.action = na.pass, plot = FALSE)$acf[2], silent = TRUE)
rho <- if (inherits(acf1, "try-error") || is.na(acf1)) 0 else max(0, min(0.9, as.numeric(acf1)))
# 2.3 Definir bloques AR.start (TRUE si hay salto grande entre meses)
mensual_total <- mensual_total |>
arrange(mes) |>
mutate(ARstart = c(TRUE, diff(mes) > 45)) # nuevo bloque si gap > 45 días
# 2.4 BAM con AR(1)
gam_total_ar <- bam(
conteo_total ~ s(t_meses, k = 20) + s(mes_num, bs = "cc", k = 12),
data = mensual_total,
family = nb(), method = "fREML",
knots = list(mes_num = c(0.5, 12.5)),
AR.start = ARstart, rho = rho
)
summary(gam_total_ar)
##
## Family: Negative Binomial(5.711)
## Link function: log
##
## Formula:
## conteo_total ~ s(t_meses, k = 20) + s(mes_num, bs = "cc", k = 12)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.44021 0.04536 97.88 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(t_meses) 4.159 5.18 2.82 0.018749 *
## s(mes_num) 4.083 10.00 2.18 0.000236 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.286 Deviance explained = 31.5%
## fREML = 138.03 Scale est. = 1 n = 91
# Serie + ajuste (AR)
df_pred_total <- pred_ci(gam_total_ar, mensual_total)
p_total <- ggplot() +
geom_line(data = mensual_total, aes(mes, conteo_total),
linewidth = 0.6, color = "grey55") +
geom_point(data = mensual_total, aes(mes, conteo_total),
size = 1.1, color = "grey35") +
geom_ribbon(data = df_pred_total, aes(mes, ymin = lwr, ymax = upr),
alpha = 0.18, fill = pal_okabe[7]) +
geom_line(data = df_pred_total, aes(mes, mu),
linewidth = 1.2, color = pal_okabe[7]) +
scale_x_date(date_breaks = "1 year", date_labels = "%Y",
expand = expansion(mult = c(.01,.03))) +
labs(
title = "Conteo total mensual y tendencia ajustada (GAM NegBin con AR(1))",
subtitle = "Línea y banda: predicción e IC95% del modelo",
x = "Mes", y = "Individuos (total mensual)"
)
p_total
ggsave_Ara("figs/01_total_ar.png", p_total)
# Efecto estacional (cíclico) — factor relativo exp(s(mes))
s_m <- smooth_curve(gam_total_ar, "s(mes_num)",
x_candidates = c("mes_num","x",".x"), out_x = "mes_num",
level = 0.95, linkinv = exp)
p_est <- ggplot(s_m, aes(mes_num, mu, group = 1)) +
geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.18, fill = pal_okabe[3]) +
geom_line(linewidth = 1.2, color = pal_okabe[3]) +
geom_point(size = 1.2, color = pal_okabe[3]) +
scale_x_continuous(breaks = 1:12, labels = c("E","F","M","A","M","J","J","A","S","O","N","D")) +
labs(title = "Efecto estacional (exp(s(mes)))", x = "Mes", y = "Factor relativo",
caption = "IC puntuales 95%")
p_est
ggsave_Ara("figs/02_estacionalidad.png", p_est)
# Derivadas simultáneas (95%) sobre la misma malla temporal
grid <- mensual_total %>% arrange(mes) %>% select(mes, t_meses, mes_num)
der_ar <- gratia::derivatives(
gam_total_ar, term = "s(t_meses)",
newdata = grid, interval = "simultaneous"
)
# Armoniza nombres posibles
if (".derivative" %in% names(der_ar)) der_ar <- der_ar %>% rename(derivative = .derivative)
if (".lower_ci" %in% names(der_ar)) der_ar <- der_ar %>% rename(lower_ci = .lower_ci)
if (".upper_ci" %in% names(der_ar)) der_ar <- der_ar %>% rename(upper_ci = .upper_ci)
der_ar <- der_ar %>% mutate(mes = grid$mes,
sig_inc = lower_ci > 0,
sig_dec = upper_ci < 0)
inc_tramos_ar <- segmentar_tramos(der_ar, "sig_inc")
dec_tramos_ar <- segmentar_tramos(der_ar, "sig_dec")
list(
n_tramos_inc = nrow(inc_tramos_ar),
n_tramos_dec = nrow(dec_tramos_ar)
)
## $n_tramos_inc
## [1] 0
##
## $n_tramos_dec
## [1] 0
4a) Sombrado de tramos significativos (95% simultáneo, AR(1)) en la figura del total
# 1) Filtra tramos "sostenidos" (puedes ajustar este mínimo)
min_meses <- 3
inc_rects <- inc_tramos_ar %>%
dplyr::filter(n_meses >= min_meses) %>%
dplyr::transmute(xmin = inicio, xmax = fin, ymin = -Inf, ymax = Inf, tipo = "Incremento")
dec_rects <- dec_tramos_ar %>%
dplyr::filter(n_meses >= min_meses) %>%
dplyr::transmute(xmin = inicio, xmax = fin, ymin = -Inf, ymax = Inf, tipo = "Decremento")
rects <- dplyr::bind_rows(inc_rects, dec_rects)
# 2) Agrega el sombreado a la figura del total (si no hay tramos, no pasa nada)
p_total_sombreado <- p_total +
(if (nrow(rects) > 0) {
ggplot2::geom_rect(
data = rects,
aes(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax, fill = tipo),
inherit.aes = FALSE, alpha = 0.12
)
} else NULL) +
scale_fill_manual(values = c("Incremento" = "#009E73", "Decremento" = "#D55E00")) +
guides(fill = guide_legend(title = NULL, override.aes = list(alpha = 0.35))) +
labs(subtitle = paste0(
"Zonas sombreadas: tramos con pendiente significativa (",
min_meses, "+ meses, 95% simultáneo, AR(1))"
))
p_total_sombreado
ggsave_Ara("figs/01_total_ar_sombreado.png", p_total_sombreado)
gam_adultos <- gam(adultos_total ~ s(t_meses, k = 20) + s(mes_num, bs = "cc", k = 12),
data = mensual_total, family = nb(), method = "REML",
knots = list(mes_num = c(0.5, 12.5)))
gam_juveniles <- gam(juveniles_total ~ s(t_meses, k = 20) + s(mes_num, bs = "cc", k = 12),
data = mensual_total, family = nb(), method = "REML",
knots = list(mes_num = c(0.5, 12.5)))
pred_total <- pred_ci(gam_total_ar, mensual_total) %>% mutate(serie = "Total")
pred_ad <- pred_ci(gam_adultos, mensual_total) %>% mutate(serie = "Adultos")
pred_ju <- pred_ci(gam_juveniles, mensual_total) %>% mutate(serie = "Juveniles")
pred_all <- bind_rows(pred_total, pred_ad, pred_ju)
obs_long <- mensual_total |>
select(mes, Total = conteo_total, Adultos = adultos_total, Juveniles = juveniles_total) |>
pivot_longer(-mes, names_to = "serie", values_to = "y")
colores <- c("Total" = pal_okabe[7], "Adultos" = pal_okabe[9], "Juveniles" = pal_okabe[6])
p_comb <- ggplot() +
geom_line(data = obs_long, aes(mes, y, color = serie), linewidth = 0.3, alpha = 0.55) +
geom_point(data = obs_long, aes(mes, y, color = serie), size = 0.7, alpha = 0.55) +
geom_ribbon(data = pred_all, aes(mes, ymin = lwr, ymax = upr, fill = serie),
alpha = 0.12, color = NA) +
geom_line(data = pred_all, aes(mes, mu, color = serie), linewidth = 0.6) +
scale_color_manual(values = colores, name = "Serie") +
scale_fill_manual(values = colores, guide = "none") +
scale_x_date(date_breaks = "1 year", date_labels = "%Y",
expand = expansion(mult = c(.01, .03))) +
labs(
title = "Guacamayas: total, adultos y juveniles",
subtitle = "Curvas GAM con IC95% (bandas) sobre series observadas",
x = "Mes", y = "Individuos por mes"
)
p_comb
ggsave_Ara("figs/03_combinado_total_adultos_juveniles.png", p_comb)
# Misma figura pero en escala relativa
# Máximo observado por serie para escalar 0–1 (evita división por 0)
max_by <- obs_long |>
dplyr::group_by(serie) |>
dplyr::summarise(mmax = max(y, na.rm = TRUE), .groups = "drop") |>
dplyr::mutate(mmax = pmax(mmax, 1))
obs_long_norm <- obs_long |>
dplyr::left_join(max_by, by = "serie") |>
dplyr::mutate(y_norm = y / mmax)
pred_all_norm <- pred_all |>
dplyr::left_join(max_by, by = "serie") |>
dplyr::mutate(mu_norm = mu / mmax,
lwr_norm = lwr / mmax,
upr_norm = upr / mmax)
p_comb_norm <- ggplot() +
geom_line(data = obs_long_norm, aes(mes, y_norm, color = serie),
linewidth = 0.3, alpha = 0.55) +
geom_point(data = obs_long_norm, aes(mes, y_norm, color = serie),
size = 0.7, alpha = 0.55) +
geom_ribbon(data = pred_all_norm, aes(mes, ymin = lwr_norm, ymax = upr_norm, fill = serie),
alpha = 0.12, color = NA) +
geom_line(data = pred_all_norm, aes(mes, mu_norm, color = serie),
linewidth = 0.6) +
scale_color_manual(values = colores, name = "Serie") +
scale_fill_manual(values = colores, guide = "none") +
scale_x_date(date_breaks = "1 year", date_labels = "%Y",
expand = expansion(mult = c(.01, .03))) +
scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
labs(
title = "Guacamayas: total, adultos y juveniles (escala relativa)",
subtitle = "Cada serie normalizada por su máximo observado (0–1)",
x = "Mes", y = "Escala relativa"
)
p_comb_norm
ggsave_Ara("figs/03a_combinado_total_adultos_juveniles.png", p_comb_norm)
p_comb_norm_noraw<-ggplot() +geom_ribbon(data = pred_all_norm, aes(mes, ymin = lwr_norm, ymax = upr_norm, fill = serie),
alpha = 0.12, color = NA) +
geom_line(data = pred_all_norm, aes(mes, mu_norm, color = serie),
linewidth = 0.4) +
scale_color_manual(values = colores, name = "Serie") +
scale_fill_manual(values = colores, guide = "none") +
scale_x_date(date_breaks = "1 year", date_labels = "%Y",
expand = expansion(mult = c(.01, .03))) +
scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
labs(
title = "Guacamayas: total, adultos y juveniles (escala relativa)",
subtitle = "Cada serie normalizada por su máximo observado (0–1)/sin valores raw",
x = "Mes", y = "Escala relativa"
)
p_comb_norm_noraw
ggsave_Ara("figs/03b_combinado_total_adultos_juveniles.png", p_comb_norm_noraw)
mensual_tr <- mensual_tr %>% mutate(transecto = forcats::fct_infreq(transecto))
gam_tr <- gam(
conteo ~ s(t_meses, by = transecto, bs = "fs", k = 10) + s(mes_num, bs = "cc", k = 12),
data = mensual_tr, family = nb(), method = "REML",
knots = list(mes_num = c(0.5, 12.5))
)
df_pred_tr <- pred_ci(gam_tr, mensual_tr)
p_tr <- ggplot() +
geom_line(data = mensual_tr, aes(mes, conteo), linewidth = 0.5, color = "grey60") +
geom_ribbon(data = df_pred_tr, aes(mes, ymin = lwr, ymax = upr), alpha = 0.16, fill = pal_okabe[3]) +
geom_line(data = df_pred_tr, aes(mes, mu), linewidth = 0.9, color = pal_okabe[1]) +
facet_wrap(~ transecto, scales = "free_y") +
scale_x_date(date_breaks = "1 year", date_labels = "%Y") +
labs(title = "Tendencias por transecto (GAM)", subtitle = "Ajuste con IC95%",
x = "Mes", y = "Individuos por mes y transecto")
p_tr
ggsave_Ara("figs/04_transectos_tendencias.png", p_tr, width = 12, height = 7)
# Curva del smooth temporal (exp(s(t))) para amplitud relativa
s_t <- smooth_curve(gam_total_ar, "s(t_meses)",
x_candidates = c("t_meses","x",".x"), out_x = "t_meses",
level = 0.95, linkinv = exp)
amplitud_pct <- 100 * (max(s_t$mu) / min(s_t$mu) - 1)
i_max <- which.max(s_t$mu); i_min <- which.min(s_t$mu)
t_peak <- s_t$t_meses[i_max]; t_trough <- s_t$t_meses[i_min]
t0 <- min(mensual_total$mes, na.rm = TRUE)
fecha_peak <- t0 %m+% months(round(t_peak))
fecha_trough <- t0 %m+% months(round(t_trough))
list(
amplitud_pct = amplitud_pct,
fecha_valle = fecha_trough,
fecha_pico = fecha_peak
)
## $amplitud_pct
## [1] 68.44102
##
## $fecha_valle
## [1] "2017-02-01"
##
## $fecha_pico
## [1] "2021-03-01"
# Construir eje de fechas para el smooth (t_meses es "meses desde t0")
s_t_plot <- s_t %>%
mutate(
fecha = t0 %m+% months(round(t_meses)) # fecha aproximada para cada t_meses
)
# Punto de pico y valle (para anotarlos)
pt_peak <- s_t_plot %>% filter(fecha == fecha_peak) %>% slice_tail(n=1)
pt_trough <- s_t_plot %>% filter(fecha == fecha_trough) %>% slice_tail(n=1)
# Gráfica
p_amplitud <- ggplot(s_t_plot, aes(x = fecha, y = mu)) +
geom_ribbon(aes(ymin = lwr, ymax = upr),
alpha = 0.20, fill = if (exists("pal_okabe")) pal_okabe[1] else "steelblue") +
geom_line(linewidth = 1.2,
color = if (exists("pal_okabe")) pal_okabe[1] else "steelblue") +
geom_hline(yintercept = 1, linetype = "dashed", color = "grey40") +
# Marcas de pico y valle
{ if (nrow(pt_peak)>0) geom_vline(xintercept = pt_peak$fecha, linetype = "dashed", color = "red") } +
{ if (nrow(pt_trough)>0) geom_vline(xintercept = pt_trough$fecha, linetype = "dashed", color = "red3") } +
{ if (nrow(pt_peak)>0) geom_point(data = pt_peak, color = "red", size = 2) } +
{ if (nrow(pt_trough)>0) geom_point(data = pt_trough, color = "red3", size = 2) } +
scale_x_date(date_breaks = "1 year", date_labels = "%Y",
expand = expansion(mult = c(.01,.03))) +
labs(
title = "Efecto multi-anual: exp(s(t_meses))",
subtitle = paste0(
"Amplitud relativa ≈ ", round(amplitud_pct, 1), "% · ",
"Pico: ", format(fecha_peak, "%Y-%m"), " · ",
"Valle: ", format(fecha_trough, "%Y-%m"),
" (línea punteada horizontal = factor 1)"
),
x = "Fecha",
y = "Factor multiplicativo (exp(s(t)))"
) +
theme_minimal()
p_amplitud
ggsave_Ara("figs/02b_Estacionalidad_multianual.png", p_amplitud)
# Corre solo si usar_lluvia_ar == TRUE
if (isTRUE(usar_lluvia_ar)) {
# --- 8.1 Lluvia mensual ÚNICA (idéntica en todos los transectos) ---
lluvia_mensual <- mensual_tr %>%
group_by(mes) %>%
summarise(
lluvia = {
v <- unique(lluvia[!is.na(lluvia)])
if (length(v) == 0) NA_real_
else if (length(v) == 1) v
else {
warning(paste(
"Mes", format(mes[1], "%Y-%m"),
"tiene múltiples valores de lluvia:",
paste(signif(v, 4), collapse = ", "),
"-> usaré el PROMEDIO para continuar."
))
mean(v, na.rm = TRUE)
}
},
.groups = "drop"
)
# --- 8.1b Agregación mensual total y unión con lluvia única ---
mensual_total <- mensual_tr %>%
group_by(mes) %>%
summarise(
conteo_total = sum(conteo, na.rm = TRUE),
adultos_total = sum(adultos, na.rm = TRUE),
juveniles_total = sum(juveniles, na.rm = TRUE),
.groups = "drop"
) %>%
left_join(lluvia_mensual, by = "mes") %>%
mutate(
mes_num = lubridate::month(mes),
t_meses = as.numeric(difftime(mes, min(mes, na.rm = TRUE), units = "days"))/30.4375
) %>%
arrange(mes) %>%
mutate(ARstart = c(TRUE, diff(mes) > 45)) # nuevo bloque si gap > 45 días
# --- 8.2 Lags de precipitación (mes actual y 2 previos) ---
add_lags <- function(df, var = "lluvia", k = 2) {
df <- df %>% arrange(mes)
for (i in 0:k) df[[paste0(var, "_lag", i)]] <- dplyr::lag(df[[var]], i)
df
}
mensual_total_ll <- mensual_total %>% add_lags("lluvia", 2)
# --- 8.3 Estimar rho: GAM sin AR (con lluvia) para ACF(1) inicial ---
gam_total_ll0 <- mgcv::gam(
conteo_total ~ s(t_meses, k = 20) + s(mes_num, bs = "cc", k = 12) +
s(lluvia_lag0, k = 5) + s(lluvia_lag1, k = 5) + s(lluvia_lag2, k = 5),
data = mensual_total_ll,
family = nb(), method = "REML",
select = TRUE,
knots = list(mes_num = c(0.5, 12.5)),
na.action = na.exclude
)
acf1_ll <- try(acf(residuals(gam_total_ll0, type = "pearson"),
na.action = na.pass, plot = FALSE)$acf[2], silent = TRUE)
rho_ll <- if (inherits(acf1_ll, "try-error") || is.na(acf1_ll)) 0 else
max(0, min(0.9, as.numeric(acf1_ll)))
# --- 8.4 Modelo con AR(1) + precipitación (rezagos 0–2) ---
gam_total_lluvia_ar <- mgcv::bam(
conteo_total ~ s(t_meses, k = 20) + s(mes_num, bs = "cc", k = 12) +
s(lluvia_lag0, k = 5) + s(lluvia_lag1, k = 5) + s(lluvia_lag2, k = 5),
data = mensual_total_ll,
family = nb(), method = "fREML",
select = TRUE,
knots = list(mes_num = c(0.5, 12.5)),
AR.start = ARstart, rho = rho_ll,
na.action = na.exclude
)
# --- 8.4bis Predicción para figura (relleno SOLO visual en lags NA) ---
typ_ll <- mensual_total_ll %>% summarise(
lluvia_lag0 = mean(lluvia_lag0, na.rm = TRUE),
lluvia_lag1 = mean(lluvia_lag1, na.rm = TRUE),
lluvia_lag2 = mean(lluvia_lag2, na.rm = TRUE)
)
newdata_ll_full <- mensual_total_ll %>%
mutate(
lluvia_lag0 = ifelse(is.na(lluvia_lag0), typ_ll$lluvia_lag0, lluvia_lag0),
lluvia_lag1 = ifelse(is.na(lluvia_lag1), typ_ll$lluvia_lag1, lluvia_lag1),
lluvia_lag2 = ifelse(is.na(lluvia_lag2), typ_ll$lluvia_lag2, lluvia_lag2)
)
pred_lluvia_ar <- pred_ci(gam_total_lluvia_ar, newdata_ll_full)
p_total_lluvia_ar <- ggplot() +
geom_line(data = mensual_total_ll, aes(mes, conteo_total),
linewidth = 0.6, color = "grey55") +
geom_point(data = mensual_total_ll, aes(mes, conteo_total),
size = 1.1, color = "grey35") +
geom_ribbon(data = pred_lluvia_ar, aes(mes, ymin = lwr, ymax = upr),
alpha = 0.18, fill = pal_okabe[2]) +
geom_line(data = pred_lluvia_ar, aes(mes, mu),
linewidth = 1.2, color = pal_okabe[2]) +
scale_x_date(date_breaks = "1 year", date_labels = "%Y",
expand = expansion(mult = c(.01,.03))) +
labs(
title = "Total con lluvia (GAM NegBin + AR(1))",
subtitle = "Línea: ajuste con AR(1); lags rellenados SOLO para visualización continua",
x = "Mes", y = "Individuos/mes"
)
p_total_lluvia_ar
ggsave_Ara("figs/01c_total_lluvia_ar.png", p_total_lluvia_ar)
# --- 8.5 Comparación con el modelo principal (AR sin lluvia) ---
aic_ar <- tryCatch(AIC(gam_total_ar), error = function(e) NA_real_)
aic_ll <- tryCatch(AIC(gam_total_lluvia_ar), error = function(e) NA_real_)
crit_ar <- suppressWarnings(gam_total_ar$gcv.ubre)
crit_ll <- suppressWarnings(gam_total_lluvia_ar$gcv.ubre)
delta_aic <- aic_ar - aic_ll # > 2 favorece el modelo con lluvia
delta_freml <- crit_ar - crit_ll # > 0 favorece el modelo con lluvia
# Significancia de los smooths de lluvia
stab_rain <- summary(gam_total_lluvia_ar)$s.table
rn <- rownames(stab_rain)
id_rain <- grepl("lluvia", rn, ignore.case = TRUE)
n_sig_rain <- if (any(id_rain)) sum(stab_rain[id_rain, "p-value"] < 0.05, na.rm = TRUE) else 0
lluvia_diag <- list(
AIC = c(sin_lluvia = aic_ar, con_lluvia = aic_ll, delta = delta_aic),
fREML = c(sin_lluvia = crit_ar, con_lluvia = crit_ll, delta = delta_freml),
smooths_lluvia_signif = n_sig_rain
)
print(lluvia_diag)
# Flags para interpretación
lluvia_mejora <- (is.finite(delta_aic) && delta_aic > 2) ||
(is.finite(delta_freml) && delta_freml > 0)
lluvia_signif <- n_sig_rain > 0
} else {
# Placeholders si no se corre lluvia
lluvia_mejora <- FALSE
lluvia_signif <- FALSE
delta_aic <- NA_real_; delta_freml <- NA_real_
}
## $AIC
## sin_lluvia con_lluvia delta
## 927.1563 194.6377 732.5186
##
## $fREML
## sin_lluvia.fREML con_lluvia.fREML delta.fREML
## 138.03029 36.59609 101.43420
##
## $smooths_lluvia_signif
## [1] 0
Nota: el if (isTRUE(usuar_lluvia_ar)) con el typo es
deliberado para que, si copias este bloque aislado, notes que debes usar
la variable usar_lluvia_ar correcta.
Puedes borrarlo si prefieres.
s_tab <- summary(gam_total_ar)$s.table
rownames(s_tab) <- make.names(rownames(s_tab))
p_t <- as.numeric(s_tab["s.t_meses.", "p-value"])
edf_t<- round(as.numeric(s_tab["s.t_meses.", "edf"]), 2)
p_m <- as.numeric(s_tab["s.mes_num.", "p-value"])
edf_m<- round(as.numeric(s_tab["s.mes_num.", "edf"]), 2)
devexp <- round(summary(gam_total_ar)$dev.expl * 100, 1)
r2adj <- round(summary(gam_total_ar)$r.sq * 100, 1)
# Derivadas simultáneas AR(1)
grid <- mensual_total %>% arrange(mes) %>% select(mes, t_meses, mes_num)
der_ar <- gratia::derivatives(gam_total_ar, term = "s(t_meses)",
newdata = grid, interval = "simultaneous")
## Use of the `newdata` argument is deprecated.
## Instead, use the data argument `data`.
if (".lower_ci" %in% names(der_ar)) names(der_ar)[names(der_ar)==".lower_ci"] <- "lower_ci"
if (".upper_ci" %in% names(der_ar)) names(der_ar)[names(der_ar)==".upper_ci"] <- "upper_ci"
n_inc <- sum(der_ar$lower_ci > 0, na.rm = TRUE) # conteo de puntos (no tramos)
n_dec <- sum(der_ar$upper_ci < 0, na.rm = TRUE)
# Amplitud multi-anual (exp(s(t)))
s_t <- smooth_curve(gam_total_ar, "s(t_meses)",
x_candidates = c("t_meses","x",".x"), out_x = "t_meses",
level = 0.95, linkinv = exp)
amplitud_pct <- 100 * (max(s_t$mu) / min(s_t$mu) - 1)
i_max <- which.max(s_t$mu); i_min <- which.min(s_t$mu)
t0 <- min(mensual_total$mes, na.rm = TRUE)
fecha_peak <- t0 %m+% months(round(s_t$t_meses[i_max]))
fecha_valle <- t0 %m+% months(round(s_t$t_meses[i_min]))
cat("### Interpretación\n\n")
cat("- **Modelo principal:** GAM NegBin con corrección **AR(1)**.\n")
cat("- **Ajuste global:** desviancia explicada **", devexp, "%**; R²(adj) **", r2adj, "%**.\n", sep = "")
cat("- **Estacionalidad anual** `s(mes)`: **p = ", format.pval(p_m), "**; edf ≈ ", edf_m, ".\n", sep = "")
s(mes): p
= 0.00023569; edf ≈ 4.08.cat("- **Tendencia multi-anual** `s(t)`: **p = ", format.pval(p_t), "**; edf ≈ ", edf_t,
" ⇒ estructura **no lineal** (ondulante) más allá de la estacionalidad.\n", sep = "")
s(t): p
= 0.018749; edf ≈ 4.16 ⇒ estructura no lineal
(ondulante) más allá de la estacionalidad.cat("- **Amplitud multi-anual** (controlando estacionalidad): ~", round(amplitud_pct, 1),
"% (valle ≈ ", format(fecha_valle, "%Y-%m"), " → pico ≈ ", format(fecha_peak, "%Y-%m"), ").\n", sep = "")
# Cambios locales robustos (derivadas simultáneas)
if (n_inc + n_dec == 0) {
cat("- **Cambios locales (derivadas simultáneas 95%)**: no se detectaron pendientes >0 o <0 robustas ⇒ ",
"oscilaciones sin tramos direccionales sostenidos.\n", sep = "")
} else {
cat("- **Cambios locales (derivadas simultáneas 95%)**: se detectaron puntos con pendiente ≠0; ",
"ver figura sombreada si segmentas tramos.\n", sep = "")
}
# Lluvia con AR(1) (si se corrió)
if (exists("usar_lluvia_ar") && isTRUE(usar_lluvia_ar)) {
if (isTRUE(lluvia_mejora)) {
cat("- **Lluvia (AR(1))**: mejora del ajuste (ΔAIC > 2 y/o ΔfREML > 0). ")
if (isTRUE(lluvia_signif)) {
cat("Al menos un término `s(lluvia_lag*)` resultó **significativo**. ",
"La señal de lluvia se considera **relevante** bajo AR(1).\n", sep = "")
} else {
cat("No obstante, los `s(lluvia_lag*)` no fueron globalmente significativos.\n", sep = "")
}
} else {
cat("- **Lluvia (AR(1))**: **no** mejoró el ajuste (ΔAIC ≤ 2 y/o ΔfREML ≤ 0) ",
"y/o los términos `s(lluvia_lag*)` no fueron significativos ⇒ ",
"**sin efecto de lluvia detectable** en esta escala.\n", sep = "")
}
}
s(lluvia_lag*) no fueron
globalmente significativos.cat("\n**Conclusión:** La población muestra **estacionalidad anual** clara y una **variación multi-anual significativa** (ondulante, no monotónica). ",
"Bajo un criterio conservador (AR(1) + bandas simultáneas), no se evidencian tramos sostenidos de aumento/disminución. ",
"La lluvia, cuando se evalúa con AR(1), ",
if (exists('usar_lluvia_ar') && isTRUE(usar_lluvia_ar) && isTRUE(lluvia_mejora)) "puede aportar algo" else "no aporta mejora apreciable",
".\n", sep = "")
Conclusión: La población muestra estacionalidad anual clara y una variación multi-anual significativa (ondulante, no monotónica). Bajo un criterio conservador (AR(1) + bandas simultáneas), no se evidencian tramos sostenidos de aumento/disminución. La lluvia, cuando se evalúa con AR(1), puede aportar algo. ## Detalles útiles
bam):
valores menores son mejores; en el código uso ΔfREML =
fREML(sin) − fREML(con) > 0 para indicar mejora.select=TRUE deja a mgcv
penalizar a cero smooths innecesarios: si la lluvia no
aporta, verás edf muy pequeñas o p grandes.library(dplyr)
library(lubridate)
library(mgcv)
library(ggplot2)
# 1) Unir pollos anual
pollos_anual <- tibble::tibble(
anio = 2013:2024,
n_pollos = c(0, 0, 0, 0, 0, 15, 22, 25,27,12,12,9) # <- tus valores reales
)
mensual_pollos <- mensual_total %>%
mutate(anio = year(mes)) %>%
left_join(pollos_anual, by = "anio")
# (Opcional) Poner 0 si faltan datos de algunos años
# mensual_pollos <- mensual_pollos %>% mutate(n_pollos = coalesce(n_pollos, 0))
# 2) Asegurar ARstart
if (!"ARstart" %in% names(mensual_pollos)) {
mensual_pollos <- mensual_pollos %>%
arrange(mes) %>%
mutate(ARstart = c(TRUE, diff(mes) > 45))
}
# 3) Filtrar filas completas
model_df <- mensual_pollos %>%
filter(!is.na(conteo_total), !is.na(n_pollos), !is.na(t_meses), !is.na(mes_num))
cat("Filas disponibles:", nrow(model_df), "\n")
## Filas disponibles: 91
# 4) Modelo base (SIN pollos), SIN lluvia
gam_base_sin_lluvia <- mgcv::bam(
conteo_total ~ s(t_meses, k = 10) + s(mes_num, bs = "cc", k = 12),
data = model_df,
family = nb(), method = "fREML",
select = TRUE,
knots = list(mes_num = c(0.5, 12.5)),
AR.start = model_df$ARstart,
na.action = na.exclude
)
# 5) Modelo con pollos (sin lluvia), efecto NO lineal para extraer curva parcial
gam_pollos <- mgcv::bam(
conteo_total ~ s(t_meses, k = 10) + s(mes_num, bs = "cc", k = 12) +
s(n_pollos, k = 5),
data = model_df,
family = nb(), method = "fREML",
select = TRUE,
knots = list(mes_num = c(0.5, 12.5)),
AR.start = model_df$ARstart,
na.action = na.exclude
)
# 6) Comparación y evidencia
aic_base <- tryCatch(AIC(gam_base_sin_lluvia), error = function(e) NA_real_)
aic_pol <- tryCatch(AIC(gam_pollos), error = function(e) NA_real_)
crit_base <- suppressWarnings(gam_base_sin_lluvia$gcv.ubre)
crit_pol <- suppressWarnings(gam_pollos$gcv.ubre)
delta_aic <- aic_base - aic_pol # > 2 favorece incluir pollos
delta_freml <- crit_base - crit_pol # > 0 favorece incluir pollos
stab <- summary(gam_pollos)$s.table
p_pollos <- if ("s(n_pollos)" %in% rownames(stab)) stab["s(n_pollos)", "p-value"] else NA_real_
cat("ΔAIC =", round(delta_aic,3),
"| ΔfREML =", round(delta_freml,3),
"| p-valor s(n_pollos) =", signif(p_pollos,3), "\n")
## ΔAIC = 4.4 | ΔfREML = 0.66 | p-valor s(n_pollos) = 0.0154
# 7) Predicción serie temporal con IC 95% (modelo con pollos)
pred_link <- predict(gam_pollos, newdata = model_df, type = "link", se.fit = TRUE)
z <- qnorm(0.975)
fit_link <- pred_link$fit
se_link <- pred_link$se.fit
lwr_link <- fit_link - z * se_link
upr_link <- fit_link + z * se_link
model_df <- model_df %>%
mutate(
fitted = gam_pollos$family$linkinv(fit_link),
lwr = gam_pollos$family$linkinv(lwr_link),
upr = gam_pollos$family$linkinv(upr_link)
)
p_ts <- ggplot(model_df, aes(x = mes)) +
geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.20) +
geom_line(aes(y = fitted), linewidth = 1.2) +
geom_point(aes(y = conteo_total), color = "grey30", size = 1.1) +
geom_line(aes(y = conteo_total), color = "grey55", linewidth = 0.6) +
scale_x_date(date_breaks = "1 year", date_labels = "%Y",
expand = expansion(mult = c(.01,.03))) +
labs(
title = "Tendencia poblacional (sin lluvia) con efecto de pollos rescatados",
subtitle = "GAM NegBin + AR(1) implícito vía bloques; IC 95% sobre la predicción",
x = "Mes", y = "Individuos/mes",
caption = paste0("ΔAIC=", round(delta_aic,2),
" | ΔfREML=", round(delta_freml,2),
" | p s(n_pollos)=", signif(p_pollos,3))
)
print(p_ts)
ggsave_Ara("figs/03_pollos_sin_lluvia_ts.png", p_ts)
# 8) Efecto parcial de 'n_pollos' con IC 95%
# Usa tu helper smooth_curve() para extraer en respuesta
idx <- which(attr(gam_pollos$smooth, "term.labels") == "s(n_pollos)")
# Crear una secuencia de valores de n_pollos para predecir el efecto
grid_pollos <- data.frame(
n_pollos = seq(min(model_df$n_pollos, na.rm=TRUE),
max(model_df$n_pollos, na.rm=TRUE),
length.out = 100),
t_meses = mean(model_df$t_meses, na.rm=TRUE),
mes_num = mean(model_df$mes_num, na.rm=TRUE)
)
# Predicciones del smooth específico
pred_smooth <- predict(gam_pollos, newdata = grid_pollos,
type = "terms", terms = "s(n_pollos)", se.fit = TRUE)
# Construir data.frame con IC 95% en escala de respuesta
z <- qnorm(0.975)
grid_pollos <- grid_pollos %>%
mutate(
fit = pred_smooth$fit[,1],
se = pred_smooth$se.fit[,1],
lwr = fit - z*se,
upr = fit + z*se
)
# Pasar a escala de respuesta
linkinv <- gam_pollos$family$linkinv
grid_pollos <- grid_pollos %>%
mutate(mu = linkinv(fit), lwr = linkinv(lwr), upr = linkinv(upr))
# Graficar efecto parcial
p_partial_manual <- ggplot(grid_pollos, aes(x = n_pollos, y = mu)) +
geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.2, fill = "#fb8c27") +
geom_line(color = "#fb8c27", linewidth = 1) +
labs(
title = "Efecto parcial de pollos rescatados",
subtitle = "GAM NegBin, IC 95%",
x = "Pollos rescatados (anual)",
y = "Efecto estimado sobre el conteo"
) +
theme_minimal()
print(p_partial_manual)
ggsave_Ara("figs/03b_tendencia_incluyendo_solo_pollos_.png", p_partial_manual, width = 7, height = 5, dpi = 300)
interpretacion <- tibble::tibble(
Indicador = c("ΔAIC (>2 favorece)", "ΔfREML (>0 favorece)", "p-valor s(n_pollos) (<0.05)"),
Valor = c(round(delta_aic, 3), round(delta_freml, 3), signif(p_pollos, 3)),
Evidencia = c(
ifelse(delta_aic > 2, "Sí", "No"),
ifelse(delta_freml > 0, "Sí", "No"),
ifelse(!is.na(p_pollos) && p_pollos < 0.05, "Sí", "No")
)
)
# Conclusión final
favor = (delta_aic > 2 || delta_freml > 0) && (!is.na(p_pollos) && p_pollos < 0.05)
interpretacion_final <- if (favor) {
"Evidencia estadística: Sí, el programa de rescate de pollos se asocia a una población más estable."
} else {
"Evidencia estadística: No hay suficiente señal clara del efecto del rescate en la población."
}
print(interpretacion)
## # A tibble: 3 × 3
## Indicador Valor Evidencia
## <chr> <dbl> <chr>
## 1 ΔAIC (>2 favorece) 4.4 Sí
## 2 ΔfREML (>0 favorece) 0.66 Sí
## 3 p-valor s(n_pollos) (<0.05) 0.0154 Sí
cat("\n", interpretacion_final, "\n")
##
## Evidencia estadística: Sí, el programa de rescate de pollos se asocia a una población más estable.
#1. Crear una secuencia de valores de n_pollos
grid_pollos <- data.frame(
n_pollos = seq(min(model_df$n_pollos, na.rm=TRUE),
max(model_df$n_pollos, na.rm=TRUE),
length.out = 200),
t_meses = mean(model_df$t_meses, na.rm=TRUE),
mes_num = mean(model_df$mes_num, na.rm=TRUE)
)
#2. Predicción del smooth específico (sólo n_pollos)
pred_smooth <- predict(gam_pollos, newdata = grid_pollos,
type = "terms", terms = "s(n_pollos)", se.fit = TRUE)
grid_pollos <- grid_pollos %>%
mutate(
fit_link = pred_smooth$fit[,1],
se = pred_smooth$se.fit[,1],
lwr = fit_link - 1.96*se,
upr = fit_link + 1.96*se
)
#3. Encontrar el primer punto donde el efecto esperado > 0
umbral_est <- min(grid_pollos$n_pollos[grid_pollos$fit_link > 0], na.rm = TRUE)
# (Opcional) Umbral conservador: cuando el límite inferior del IC > 0
umbral_conservador <- min(grid_pollos$n_pollos[grid_pollos$lwr > 0], na.rm = TRUE)
cat("Umbral estimado (efecto esperado > 0):", umbral_est, "\n")
## Umbral estimado (efecto esperado > 0): 12.34673
cat("Umbral conservador (IC 95% > 0):", umbral_conservador, "\n")
## Umbral conservador (IC 95% > 0): 12.34673
p_partial_positive<-ggplot(grid_pollos, aes(x = n_pollos, y = fit_link)) +
geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.2, fill = "#fb8c27") +
geom_line(color = "#fb8c27", size = 1) +
geom_hline(yintercept = 0, linetype = "dashed") +
geom_vline(xintercept = umbral_est, linetype = "dotted", color = "red") +
labs(x = "Pollos rescatados (anual)",
y = "Efecto parcial (escala log)",
title = "Umbral de rescate de pollos para efecto positivo")
ggsave_Ara("figs/05_Umbral_rescate_efecto_positivo.png", p_partial_positive)
10.1) Transformando a conteos esperados en relación al número de pollos rescatados
#4. Crear grid de n_pollos como antes
grid_pollos <- data.frame(
n_pollos = seq(min(model_df$n_pollos, na.rm = TRUE),
max(model_df$n_pollos, na.rm = TRUE),
length.out = 100),
t_meses = mean(model_df$t_meses, na.rm = TRUE),
mes_num = mean(model_df$mes_num, na.rm = TRUE)
)
#5. Predicción de términos y error estándar para s(n_pollos)
pred_smooth <- predict(gam_pollos, newdata = grid_pollos,
type = "terms", terms = "s(n_pollos)", se.fit = TRUE)
#6. Incorporar predicción del intercepto (constante del modelo)
pred_intercept <- predict(gam_pollos, newdata = grid_pollos,
type = "link", exclude = "s(n_pollos)")
#7. IC 95% en escala link
z <- qnorm(0.975)
grid_pollos <- grid_pollos %>%
mutate(
fit_link = pred_smooth$fit[,1] + pred_intercept,
se = pred_smooth$se.fit[,1],
lwr_link = fit_link - z*se,
upr_link = fit_link + z*se
)
#8. Pasar todo a escala de conteos reales
linkinv <- gam_pollos$family$linkinv
grid_pollos <- grid_pollos %>%
mutate(
mu = linkinv(fit_link),
lwr = linkinv(lwr_link),
upr = linkinv(upr_link)
)
#9. Graficar efecto en escala de conteos
p_partial_counts <- ggplot(grid_pollos, aes(x = n_pollos, y = mu)) +
geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.2, fill = "#fb8c27") +
geom_line(color = "#fb8c27", linewidth = 1) +
labs(
x = "Pollos rescatados (anual)",
y = "Conteo esperado de individuos/mes",
title = "Efecto estimado de pollos rescatados (escala de conteos reales)",
subtitle = "Curva suavizada con IC 95%"
) +
theme_minimal()
print(p_partial_counts)
ggsave_Ara("figs/05b_Umbral_rescate_efecto_positivo_conteos.png", p_partial_counts)
10.2) ¿Añadir lag de 1 año tiene sentido?
# 1) Crear lag de 1 año (12 meses)
mensual_pollos <- mensual_pollos %>%
arrange(mes) %>%
mutate(n_pollos_lag1 = lag(n_pollos, 12)) # rezago anual
# 2) Filtrar filas completas
model_df_lag1 <- mensual_pollos %>%
filter(!is.na(conteo_total), !is.na(n_pollos), !is.na(n_pollos_lag1),
!is.na(t_meses), !is.na(mes_num), !is.na(ARstart))
cat("Filas disponibles para modelo lag 0 y lag 1 año:", nrow(model_df_lag1), "\n")
## Filas disponibles para modelo lag 0 y lag 1 año: 79
# 3) Ajustar el modelo GAM
gam_pollos_lag1 <- mgcv::bam(
conteo_total ~ s(t_meses, k = 10) +
s(mes_num, bs = "cc", k = 12) +
s(n_pollos, k = 5) + # efecto actual
s(n_pollos_lag1, k = 5), # efecto con 1 año de rezago
data = model_df_lag1,
family = nb(), method = "fREML",
select = TRUE,
knots = list(mes_num = c(0.5, 12.5)),
AR.start = model_df_lag1$ARstart,
na.action = na.exclude
)
summary(gam_pollos_lag1)
##
## Family: Negative Binomial(4.426)
## Link function: log
##
## Formula:
## conteo_total ~ s(t_meses, k = 10) + s(mes_num, bs = "cc", k = 12) +
## s(n_pollos, k = 5) + s(n_pollos_lag1, k = 5)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.4313 0.0549 80.71 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(t_meses) 1.529e-07 9 0.000 0.84912
## s(mes_num) 2.595e+00 10 1.100 0.00522 **
## s(n_pollos) 7.184e-01 4 0.638 0.04514 *
## s(n_pollos_lag1) 9.362e-01 4 0.487 0.11286
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.222 Deviance explained = 22.6%
## fREML = 114.26 Scale est. = 1 n = 79
# Predicciones parciales para ambos smooths
grid_pollos0 <- data.frame(
n_pollos = seq(min(model_df_lag1$n_pollos, na.rm = TRUE),
max(model_df_lag1$n_pollos, na.rm = TRUE),
length.out = 100),
n_pollos_lag1 = mean(model_df_lag1$n_pollos_lag1, na.rm = TRUE),
t_meses = mean(model_df_lag1$t_meses, na.rm = TRUE),
mes_num = mean(model_df_lag1$mes_num, na.rm = TRUE)
)
# Grid para n_pollos_lag1: variamos lag1, fijamos n_pollos
grid_pollos1 <- data.frame(
n_pollos_lag1 = seq(min(model_df_lag1$n_pollos_lag1, na.rm = TRUE),
max(model_df_lag1$n_pollos_lag1, na.rm = TRUE),
length.out = 100),
n_pollos = mean(model_df_lag1$n_pollos, na.rm = TRUE),
t_meses = mean(model_df_lag1$t_meses, na.rm = TRUE),
mes_num = mean(model_df_lag1$mes_num, na.rm = TRUE)
)
# Predicciones para smooth de n_pollos
pred0 <- predict(gam_pollos_lag1, newdata = grid_pollos0,
type = "terms", terms = "s(n_pollos)", se.fit = TRUE)
grid_pollos0 <- grid_pollos0 %>%
mutate(fit = pred0$fit[,1], se = pred0$se.fit[,1],
lwr = fit - 1.96*se, upr = fit + 1.96*se)
# Predicciones para smooth de n_pollos_lag1
pred1 <- predict(gam_pollos_lag1, newdata = grid_pollos1,
type = "terms", terms = "s(n_pollos_lag1)", se.fit = TRUE)
grid_pollos1 <- grid_pollos1 %>%
mutate(fit = pred1$fit[,1], se = pred1$se.fit[,1],
lwr = fit - 1.96*se, upr = fit + 1.96*se)
# Solo calculamos para lag1
umbral_est_1 <- min(grid_pollos1$n_pollos_lag1[grid_pollos1$fit > 0], na.rm = TRUE)
umbral_conf_1 <- min(grid_pollos1$n_pollos_lag1[grid_pollos1$lwr > 0], na.rm = TRUE)
# Pasar a escala respuesta
p0 <- ggplot(grid_pollos0, aes(x = n_pollos, y = fit)) +
geom_ribbon(aes(ymin = lwr, ymax = upr), fill = "#fb8c27", alpha = 0.2) +
geom_line(color = "#fb8c27", linewidth = 0.8) +
geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
geom_vline(xintercept = umbral_est, color = "red", linetype = "dotted") +
geom_vline(xintercept = umbral_conservador, color = "darkred", linetype = "dotted") +
labs(
x = "Pollos rescatados (año)",
y = "Efecto parcial (log)",
title = "Efecto parcial pollos rescatados",
subtitle = paste0("Umbral estimado: ", round(umbral_est,1),
" | Umbral IC95%: ", round(umbral_conservador,1))
) +
theme_minimal()
# --- Gráfico lag 1 (escala log)
p1 <- ggplot(grid_pollos1, aes(x = n_pollos_lag1, y = fit)) +
geom_ribbon(aes(ymin = lwr, ymax = upr), fill = "#04cdc1", alpha = 0.2) +
geom_line(color = "#04cdc1", linewidth = 0.8) +
geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
geom_vline(xintercept = umbral_est_1, color = "red", linetype = "dotted") +
geom_vline(xintercept = umbral_conf_1, color = "darkred", linetype = "dotted") +
labs(
x = "Pollos rescatados (año anterior)",
title = "Efecto parcial pollos lag 1 año",
subtitle = paste0("Umbral estimado: ", round(umbral_est_1,1),
" | Umbral IC95%: ", round(umbral_conf_1,1))
) +
theme_minimal()
print(p0)
print(p1)
library(patchwork)
combo <- p0 + p1 + plot_layout(widths = c(1,1)) + plot_annotation(
title = "Efecto parcial de pollos rescatados",
subtitle = "Año de rescate vs. Año anterior (IC 95%)"
)
combo
ggsave_Ara("figs/05c_Umbral_rescate_efecto_positivo_lags.png", combo)